The force network ensemble for granular packings 
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For packings of hard but not perfectly rigid particles, the length scales that govern the packing 
geometry and the contact forces are well separated. This separation of length scales is explored in the 
force network ensemble, where one studies the space of allowed force configurations for a given, frozen 
contact geometry. Here we review results of this approach, which yields nontrivial predictions for 
the effect of packing dimension and anisotropy on the contact force distribution P(f), the response 
to overall shear and point forcing, all of which can be studied in great numerical detail. Moreover, 
there are emerging analytical approaches that very effectively capture, for example, the form of force 
distributions. 

PACS numbers: 45.70.Cc, 05.40.a, 46.65.+g 



Force networks are a striking feature of granular me- 
dia [IH5] — see Fig. [I] This organization of the contact 
forces between individual grains has fascinated physicists 
for decades [6}{9]. Why are these contact forces interest- 
ing? First, the fluctuations of the forces appear to be 
unexpectedly strong, with early measurements [lOj [TT] 
and models [12] indicating that the probability distribu- 
tion function of the contact forces, is not narrowly 
distributed about its mean, but instead decays exponen- 
tially at large forces. Second, the spatial organization 
of the strong contact forces in so-called force networks 
plays an important role in the memory effects exhibited 
by granular media [13] [14]. Third, predicting the me- 
chanical properties of granular assemblies is a central goal 
of granular theory [T5H2T] , and these stresses ultimately 
originate from the contact forces and their organization. 

From a microscopic point of view, the contact forces 
arise from deformations of the constituent grains. This 
perspective is useful when developing a detailed numer- 
ical recipe, but theoretically unwieldy. Because granu- 
lar media are highly disordered the details of the con- 
tact geometry play a crucial role, response to a load is 
not affine, and macroscopic elasticity cannot be inferred 
directly from a microscopic force law [5]. Additionally, 
grain deformations are typically very small — a steel ball 
of 1 mm deforms only by a nanometer under its own 
weight, and buried under a pile of beads one meter thick, 
typical deformations are of order of 100 nm. This leads 
to a strong separation of the scale governing the geom- 
etry of the packing and the scale governing the contact 
forces. 

In this paper we review the Force Network Ensemble 
(FNE) [22], which is a model that relies on the separa- 
tion of scales relevant for the contact forces (nanometers) 
and the particle scale (millimeters). The central idea is 
that, for a given packing geometry, many different micro- 
scopic configurations of noncohesive forces satisfy force 
(and possibly torque) balance on each grain, while also 
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FIG. 1: Illustration of force networks studied in the Force Net- 
work Ensemble. Lines between disks indicate contact forces; 
larger forces have thicker lines. For a single ordered (a) or 
disordered (b) contact network, many different configurations 
of noncohesive contact forces, or force networks, satisfy force 
balance on every grain. The FNE comprises all balanced force 
networks on a fixed contact network. (From [22], ) 



satisfying boundary conditions in terms of the applied 
stresses. In other words, the forces can be seen as under- 
determined. As a simple example of force indeterminacy, 
one may think of the forces acting on a ladder that rests 
against a wall under an angle — a range of contact forces 
is possible to keep the ladder in balance. This indetermi- 
nacy carries over to collections of grains, as is illustrated 
in Fig. [I] Both ordered and disordered packings allow 
for many different force networks that respect mechan- 
ical equilibrium on each grain. By averaging over all 
possible force configurations the FNE provides a model 
for statistical properties of force networks. Of course, in 
a real physical system, the actual forces are selected by 
the history and elasticity of the ladder or particles. 



Strictly speaking, this indeterministic point of view 
only makes sense in the limit of very hard contacts. A 
subtle point is then that for perfectly rigid frictionless 
spheres, the geometry completely specifies the forces — 
frictionless hard spheres in dimension d organize such 
that their mean number of contacts z reaches a well de- 
fined limit Zi SO = 2<i, termed isostatic, such that the num- 
ber of contact forces and mechanical constraints precisely 
balance [23-26 . In more modern language, the small 
relative deformations of the particles imply that granu- 
lar media are close to the "jamming point" [5j I27H29] . 
Frictionless spheres are isostatic in this limit. Frictional 
particles, however, are generally not isostatic at the jam- 
ming point [30-33 . The idea of sampling all force con- 
figurations compatible with force balance/torque balance 
and boundary conditions is therefore on firmer footing 
physically in frictional packings. Having said that, the 
ensemble is also mathematically well-defined for friction- 
less systems, and many of the examples discussed below 
are frictionless. 

An additional motivation to study the FNE is some- 
what more abstract. The idea of sampling all possible 
configurations goes back to Edwards, who advocated con- 
sidering ensembles of all grain configurations consistent 
with some set of boundary conditions [34] . In general this 
is hard or impossible to do explicitly. The FNE can be 
seen as a restricted Edwards ensemble [17 , with an ex- 
tensive measure of the stress playing a role analogous to 
energy in equilibrium statistical mechanics [35j [36] . The 
ensemble then allows one to explore the consequences of 
and limits to the Edwards approach. For example, we 
will see that one popular notion, namely that entropy 
maximization implies robustly exponential force statis- 
tics [37 , does not survive confrontation with the FNE. 

The scope of this paper is to review the literature on 
the FNE and closely related work [22j [32j [35] HE EMS)] . 
The outline of this paper is as follows. We first moti- 
vate the ensemble in more detail and quantify the degree 
of force indeterminacy in Sec. [I] Statistical properties of 
the ensemble are discussed in Sec. [TTl with a focus on the 



contact force distribution P(f). Section III addresses 
the mechanical response to applying an external load, 
such as a uniform shear stress or a localized point force. 
The paper closes with a discussion on the successes and 
limitations of the FNE, where we also point out future 
directions. Throughout the paper we confront the ensem- 
ble predictions to experimental or numerical data, when 
available. 



I. MOTIVATION: FORCE INDETERMINACY 

The FNE crucially relies on force indeterminacy, mean- 
ing that the equations of local force and torque balance 
do not uniquely determine the forces. It is therefore 
instructive to specify it further. As an example, we 
consider a rigid ball in a groove with opening angle 20 
and contact forces f = (fn\ fn\ f[ 2) ) (see Fig.B. 
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FIG. 2: (a) A rigid ball in a groove [5 1 38 , 39 can be supported 
by purely normal forces at the wall, (b) Introducing tangen- 
tial, i.e frictional, forces allows a range of balanced force con- 
figurations. (Adapted from [5].) (c) Similar rearrangements, 
called "wheel moves" [40], do not change the net vector force 
on a grain in the frictionless triangular lattice, (d) Rearrange- 
ments in disordered packings are delocalized. Here, changes 
to normal forces (color maps to sign, thickness to magnitude) 
of a rearrangement in a frictional packing; changes to tangen- 
tial forces not shown. (From [32].) 



This system illustrates many salient features of indeter- 
minacy [5j [38] [39] . By inspection, the weight of the 
ball —mgz can be supported by purely normal forces 
fo = (c, 0,c, 0), where c = ^mg/sinO. This solution is 
not unique, however; by adding normal and tangential 
forces proportional to 5f = (— sin#, cos#, — sin#, cos#), 
other solutions of the form f + 5f can be identified. 

Indeterminacy carries over to packings of many grains, 
with the frictionless triangular lattice being an instruc- 
tive example. There each grain participates in z = 6 con- 
tacts, so there are z/2 distinct contact forces per grain. 
As each grain brings d = 2 force balance constraints, we 
anticipate (z/2) — d = 1 degree of freedom per grain, 
i.e. one rearrangement of the forces that respects force 
balance, as with the ball in a wedge. This rearrangement, 
identified in Ref. [40 , is called a wheel move (Fig. [2^) . 
The idea is that by decreasing all six forces that con- 
tact a certain grain (black arrows) while simultaneously 
increasing the six forces that lie on a shell around the 
central grain by the same amount (red arrows), the total 
vector force on all grains remains invariant. Note that 
one should be careful that all contact forces remain non- 
cohesive. 

Wheel moves can be generalized for packings with dis- 
order and friction [32] [40] [46]. In this case the rear- 
rangements of the forces generally cease to be localized, 
as demonstrated in Fig. [2]i. Just as in the frictionless 
triangular lattice, the number of independent force re- 
arrangements N w in a packing of N grains is given by 
the excess of force components over force/torque balance 
constraints, 
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Here Az = z — Z[ so and df is the number of force compo- 
nents per contact. The isostatic contact number Z{ so = 2d 
(resp. d+1) and df = 1 (resp. d) in frictionless (frictional) 
sphere packings [5]. The correction in Eq. (HI depends 
on details of the boundary conditions. 

Refs. [40] and [46] give a prescription for constructing 
a set of N w rearrangements With these in hand, 

any force network f = {fij} on a given contact network 
can be expressed by giving the weight on each force 
rearrangement: 

f = f + X> fe (5f fe . (2) 

k=l 

fo is any balanced force network with the appropriate 
stress tensor a, i.e. a particular solution. The stress ten- 
sor can be expressed in terms of the microscopic forces 
via [22] 

Here V is the volume of the contact network and 
points from the center of grain i to the center of grain j. 

The force rearrangements {Jf^} alter forces in the sys- 
tem without violating mechanical equilibrium, and are 
therefore the mechanism of force fluctuations in the FNE. 
They can be employed as Monte Carlo moves to sample 
the space of force networks [40] [47] , which can be achieved 
by varying the {wk}- It is important to note that the 
weights {wk} are strongly constrained by inequalities. 
These express the constraints that (i) all normal forces 
are noncohesive, i.e. f n > 0, and (ii) all tangential forces 
respect Coulomb's constraint |/ t | < jaf n . 

As the FNE employs a flat measure, meaning all valid 
force networks are given equal statistical weight, all the 
essential physics of the FNE is encoded in the geome- 
try of a high-dimensional space. The geometry, in turn, 
is determined by force balance and the positivity and 
Coulomb constraints. One can think of each force net- 
work as occupying a point in a space with each weight 
Wk describing an axis. In this space the noncohesive 
and Coulomb constraints each describe planar bound- 
aries where a normal force is zero or a tangential force is 
fully mobilized, respectively. These boundaries enclose a 
convex poly tope, which we call Tc and which contains 
all possible force networks {f^} on a contact network C 
[32] EH [39]. 

Indeterminacy is intimately connected with the size 
and shape of the space of force networks Tc- Globally, 
this space has dimension D = N w . Locally, the "size" of 
the space in a particular direction can be quantified by 
seeking the largest and smallest values a particular con- 
tact force component can take on, by varying the {wk} 
to the extremes of Tc [39] HE]- E.g. for the tangential 
force at contact c, the indeterminacy r\ is 

( fc\max ( £ c^min 

n( fc\ = [111 [hi / 4 x 

iyjt> i [(/ c )m ax + (/ c )m i„] > y> 



The denominator serves simply to normalize by a typical 
force scale in the packing, hence other choices are pos- 
sible. McNamara and Herrmann, who study packings 
under gravity using Contact Dynamics (CD) and MD, 
use the average weight of a grain (mg) [39]. It is appar- 
ent from Fig. ^jp that force fluctuations have a rich spatial 
structure, and indeed Ref. [39 finds that 77 is broadly dis- 
tributed and is positively correlated with contacts that 
carry large forces. Note that, although there is an up- 
per bound on the maximum normal force / n (and hence 
| ft | ) in a finite size packing [40] [47] , that bound grows 
with system size. It may therefore be useful to study the 
behavior of rj in the thermodynamic limit. 

Averaging 77 of Eq. Q over all normal and tangen- 
tial contact force components gives a purely geometric 
global measure of indeterminacy, related to the volume of 
Tc- As shown in Fig. [3] several alternative measures dis- 
play qualitatively similar, nontrivial dependence on fric- 
tion coefficient \i in frictional packings generated by CD, 
which treats perfectly rigid grains [32] 08]. Indetermi- 
nacy must be zero when the system is isostatic, which 
generically occurs for \i = and \i — >• 00 [5 . For fi- 
nite \i the global indeterminacy (77) displays a maximum 
for ji « 0.1 due to a balance between two competing ef- 
fects [32 . Increasing fi opens the Coulomb cone, which 
increases the volume of Tc without changing its dimen- 
sion. In contrast, increasing \i lowers the contact number 

and hence D. For a number of numerical protocols 
z(fi) decreases abruptly around fi w 0.1. Thus the force 
indeterminacy (77) displays a sharp signature of the pack- 
ing structure (Fig. [3]). 

II. STRESS STATISTICS 

The statistics of local measures of the stress, such as 
the force / at a contact or the pressure p on a grain, pro- 
vide a fundamental microscopic characterization of the 
material's stress state. The FNE has turned out to be an 
ideal model to probe these probability distributions. 

We first provide a brief review of characteristic proper- 
ties of the stress statistics of static packings. The second 
part of this section addresses theoretical aspects of statis- 
tics in the FNE, followed by a discussion. 

A. Characteristic features of P(f) 

Let us begin with a numerical characterization of the 
force probability density P(f) from Molecular Dynamics 
(MD) simulations on frictionless systems. Figure [4] de- 
picts P(f) for a range of confining pressures n, effectively 
changing the coordination number, for both Hookean and 
Hertzian interactions [47]. Several features stand out. 
All distributions have (i) a peak near the mean force 
(/), (ii) a nonvanishing weight as / — > + , and (Hi) a 
width comparable to (/). The latter two properties re- 
flect heterogeneity in the force network, while it has been 
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FIG. 3: Dependence of the force indeterminacy r\ on friction 
coefficient /i for contact networks equilibrated under Contact 
Dynamics 48 . Several different measures of force indeter- 
minacy show the same qualitative evolution with \i. r\\ is 
related to the relative fluctuations of contact force compo- 
nents when valid networks are randomly sampled, while r\i 
depends on the mean Euclidean distance between randomly 
sampled networks. See Ref. [38] for precise definitions. 773 is 
given by Eq. Q. (inset) The peak in 77 corresponds with a 
sharp decrease in mean contact number z. (From [48].) 
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FIG. 4: P(f) for Molecular Dynamics simulations of soft 
spheres in 2D (solid curves) and the FNE on the friction- 
less triangular lattice (dashed curves). / is normalized to (/). 
MD grains interact via one-sided (a) Hookean or (b) Hertzian 
springs, (insets) MD distributions differ in their confining 
pressure II, which tunes the contact number z relative to the 
isostatic value z- lso — 4. P(f) in the FNE shows strong quali- 
tative similarity to P(f) in MD for both force laws. (Adapted 
from H7l.) 



suggested that the peak is symptomatic of the jammed 
state [51] . These three features, the "look and feel" of 
P(f), are observed in a variety of simulations and ex- 
periments pll4l[TQ l[TT| [3Q l l5TH58]. This hints at generic 
mechanisms that can be probed within the FNE, in par- 
ticular since the qualitative features of P(f) in packings 
are insensitive to the force law. 

The simplest system to which one can apply the FNE 
is the two-dimensional contact network with triangular 
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FIG. 5: P(f) in the frictionless triangular lattice subject to a 
shear stress displays a quasi-exponential intermediate regime 
that grows with increasing anisotropy r. The distribution 
approaches a pure exponential (gray curve) in the asymptotic 
limit t — y 1. (Adapted from [3D].) 



lattice symmetry and isotropic stress [22] [40] [46], as 
in Fig. [T^i. The corresponding P(f) is represented by 
the dashed curves in Fig. [4] Indeed, the FNE predic- 
tion agrees very well with the results from Molecular 
Dynamics simulations and bears the three characteris- 
tic features of P(f). In addition, very similar results 
were obtained when numerically sampling FNE for dis- 
ordered disk packings with varying average contact num- 
bers ranging from z = 6 to 4.3 [46]. Closer approach to 
the isostatic point Z[ so = 4, where the space of allowed 
force networks vanishes, is numerically impractical; hence 
the FNE is not the appropriate tool to probe P(f) near 

Ref. [40] explores the consequences of stress anisotropy 
in the triangular lattice. Anisotropy corresponds to r > 
when the principle stresses G\ > <r 2 are unequal: 



(5) 



Anisotropy introduces a separate force distribution for 
each lattice direction. P(f) for forces aligned with 
the principal stress direction broadens and develops a 
regime of exponential decay before crossing over to faster- 
than-exponential decay (Fig. [5J. The width of this 
regime grows with increasing r. For sufficiently strong 
anisotropy, P(f) averaged over lattice directions loses its 
peak near (/), suggesting the presence of a peak in P(f) 
is not a robust signature of the jammed state [35 , 40 . In 
an ordered packing in the limit r — » 1, force is restricted 
to one lattice direction, the system is quasi-one dimen- 
sional, and P(f) becomes exactly exponential [40]. A 
similar exponential regime is observed for disordered con- 
tact networks under shear, before crossing over to asymp- 
totically Gaussian decay [46] . 

Recently, numerical results for the FNE have been ob- 
tained using umbrella sampling, with which the FNE 
can be sampled far more efficiently — and hence P(f) 
can be determined far more accurately — than ensem- 
bles generated by experiments or MD [36] [46] [47] . With 
this accuracy, which permits sampling P(f) over tens of 
decades in the tail (Fig. [6^1), it has been put beyond any 
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FIG. 6: P(f) from umbrella sampling in the frictionless tri- 
angular lattice (2D), frictionless fee lattice (3D), and a disor- 
dered packing in four dimensions (4D) with contact number 
z ~ 20.9. The dashed box indicates the plotting range in 
Fig. [4] (a) logP(/) bends downward when plotted versus /, 
indicating faster than exponential decay. If P(f) ~ exp (—f b ) 
then f~ b log P(f) approaches a constant for large /. We ob- 
serve (b) b = 2.0(1) (Gaussian decay) in 2D, (c) b — 1.6(1) in 
3D, and (d) b = 1.4(1) in 4D. (Adapted from |3B].) 



doubt that P(f) in the FNE decays faster than expo- 
nentially. This is already suggested by Figs. |4j in which 
P(f) clearly bends downwards on a semi-log plot. P(f) 
in the FNE is numerically determined in Ref. [46] to de- 
cay as exp (— / 6 ), with the exponent b dependent on the 
dimension. Consistent with theoretical arguments (see 
below), b = 2 for two dimensions (Fig.^), i.e. the tail of 
P(f) is Gaussian. For higher dimensions P(f) decays as 
a compressed exponential (Fig. |6j2,d), with b ~ 1.6(1) in 
3D and b ~ 1-4(1) in 4D. Numerics show no dependence 
of b on disorder for Az as small as 0.3 [46 , 47 . In light of 
early experiments finding exponential tails, these results 
are surprising; we return to this point below. 



To summarize, P(0) > for all frictionless force net- 
works in the FNE. Though there is often a peak for finite 
/, it can be destroyed by strong anisotropy. The force 
distribution is always wide and its asymptotic tail is set 
only by the dimension. Strong anisotropy, however, can 
open an intermediate regime of exponential decay. 



B. P(f) in theory 

Two of the robust features of P(f) identified above 
are the finite value of P(f) for vanishing force, and the 
asymptotic form of its tail. We will discuss now how, in 
the FNE, the origins of the former remain an open que- 
sion, while the latter follows from entropy maximization 
in the presence of an unanticipated constraint. 

It has been suggested that, in the presence of a flat 
measure, the finite value of P(f) for vanishing force can 
be traced to the geometric properties of the high dimen- 
sional space Tc [35]. In frictionless systems, states with 
a zero contact force sit on one of the boundary facets of 
Tc, and high dimensional spaces have the curious but 
well known property that the overwhelming majority of 
the space is close to the boundary. This fact seems to 
motivate the presence of many vanishingly small forces. 
In preparing this review, however, it became apparent 
that this reasoning does not stand up to closer scrutiny. 
Simply by dimensional analysis, the fraction v(Sf) of the 
volume of Tc within Sf of the surface must be v ~ ^"ffj- 
If this volume were simply divided among the N c contact 
forces, then P(0) « v(Sf)/N c Sf would vanish in the ther- 
modynamic limit. The fact that it does not means that a 
typical force network within Sf of one boundary facet is 
also within Sf of a finite fraction of the other boundary 
facets. Precisely how this can be understood, and how, 
e.g., P(0) depends on contact number z, is an interesting 
and open question. 

We have seen that P(f) in two dimensions has a Gaus- 
sian tail (Fig. [6J3) . To understand how this comes about, 
it is instructive to study a related distribution, namely 
that of local pressures p. The grain scale "pressures" 
Pi ~ J2j(fn)ij are convenient both because the con- 
straints of force balance enter at the grain scale and be- 
cause the pressure is a slightly coarse-grained stress. The 
form of the probability distribution P(p) is motivated in 
Ref. [36 by an entropy maximization argument, the crux 
of which we sketch here. The argument relies on the 
observation that force networks in two dimensions ad- 
mit a reciprocal representation in which forces from the 
network are used to construct tiles that tessellate space 
[59] . A tile is formed by graphically summing the contact 
forces on a grain in a right-hand fashion; Fig. [7] gives an 
example. The key observation is that, because stress a 
and volume V are both fixed in the FNE, the area of the 
tiling A = (det(j)V is also invariant [33 [49j [50] . This is 
a collective effect; the tiling only exists when every grain 
is in force balance. 

For concreteness, consider a frictionless triangular lat- 
tice with isotropic stress and pressure, a := (\/3/6r)II 1, 
where r is the grain radius. Because the stress is fixed, 
the grain-scale pressures {pi} obey a sum rule: 



5>* = mv. 



(6) 



II is therefore the average local pressure. The new ob- 
servation is that because there is a tiling and its area is 
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invariant, the local areas {a^} also obey a sum rule 

5> = V5(£)V (7) 

Note that Eq. ([7]) is satisfied automatically when Eq. (|6| 
and force balance on every grain are imposed. In ana- 
lytical calculations one typically resorts to treating local 
force balance approximately, with the consequence that 
Eq. ([7]) is violated - see Ref. [50] for a detailed discussion 
of this point. In this case, imposing the area sum rule as a 
constraint reintroduces a necessary consequence of local 
force balance. The surprise is that by incorporating this 
global constraint, instead of the O(N) local force balance 
constraints, one can successfully predict stress statistics 
in the FNE via an entropy maximization calculation [36 . 

The resulting distribution of local pressures in the fric- 
tionless triangular lattice is 

P(p) = Z~ y e" ap-7<a(p)> . (8) 

The exponent in the prefactor depends on contact net- 
work; v = 3 for the for the frictionless triangular lattice. 
Z, a, and 7 are not free parameters but Lagrange multi- 
pliers determined by normalization of P(p) and Eqs. (|6| 
and ([7]). (a(p)) oc p 2 is the average area of a tile with cor- 
responding pressure p; it appears because the area sum 
rule has been imposed. For asymptotically large pres- 
sures the quadratic term dominates and the tail of P(p) 
is Gaussian. For small 7, however, the Gaussian form 
may only become apparent deep in the tail. Eq. ([8| is in 
excellent agreement with numerics; see Fig. [7]i. Though 
it has not been shown that a Gaussian tail in P(p) re- 
quires a Gaussian tail in P(f), it is an empirical fact that 
the tails of P(f) and P(p) always have the same form. 
Hence the invariant tiling area presumably also explains 
the tail of P(f). 

The pressure sum rule, Eq. (|6|, is reminiscent of en- 
ergy in a microcanonical ensemble, which also obeys a 
sum rule. This similarity has provoked a number of au- 
thors to explore parallels between equilibrium statistical 
mechanics and stress-based ensembles of packings, un- 
der the assumption that they are also entropy maximiz- 
ing, see e.g. Refs. [T8| [6QH65] . In a statistical mechanics 
framework, a natural first step is to consider the analog 
of an ideal gas, i.e. completely neglecting correlations be- 
tween particles and maximizing entropy in the presence 
of the constraint Eq. ([6| . This "ideal gas approach" pre- 
dicts exponential tails, and is therefore not adequate to 
describe statistics of forces in the FNE. Though the cal- 
culation of Ref. [36] sketched above also neglects spatial 
correlations, by incorporating the tiling area constraint 
it retains information that the ideal gas approach throws 
out. It is remarkable that including just one more global 
constraint, rather than the many local force balance con- 
straints, so thoroughly captures the numerical distribu- 
tion (Fig.[7]i). Neglecting spatial correlations is no longer 
reasonable in the presence of a diverging length scale near 



isostaticity [66] [67] , so Eq. ([8| makes no prediction for 
this case. 

Is there a counterpart to the tiling constraint for higher 
dimensions? One can indeed construct polyhedra analo- 
gous to the polygonal tiles of Fig. [7] [36 , 68 , but whether 
the sum of their volumes is conserved remains open. If 
so, one would have P(f) ~ exp (—f b ) with b = d/(d— 1). 
Encouragingly, this is in reasonable agreement with the 
results of Fig. [6] 



C. Comparison to experiment and simulation 

How well do these features conform with P(f) from ex- 
perimental or simulated systems? Finite P(0) is indeed a 
robust feature of bulk measurements. Similarly, sampled 
distributions are broad, and a qualitative change in P(f) 
has been observed in the only experiments to systemat- 
ically vary anisotropy [4 . It is less clear if any general 
statements can be made about the form of the tail of 
P(f) in real granular systems [21 El II!B III1 ESI EIB S3 ED- 

BS1E2HZ2]. 

The earliest measurements of P(f) — made at the 
boundary of packings using the carbon paper method 
[TO] ITT] - displayed unambiguously exponential tails. 
Accompanied by the q-model [T2] . which predicts expo- 
nentials, these experiments established the expectation 
that P(f) should decay exponentially. Forces in the bulk 
are difficult to access experimentally, and there have been 
few measurements [4] [55] [57] [58] . Those there are raise 
doubts regarding the robustness of exponentials; distri- 
butions from isotropic 2D photoelastic systems [4] and 3D 
emulsions [57, 58 show downward curvature on a semi- 
log plot, suggesting faster than exponential decay. Statis- 
tics are limited, however, and it is difficult to determine 
the asymptotic form of the tail. 

Simulation results are inconsistent [2 , 29] [30] [47] I5T1 
154] [69] [70] [73] . Though some distributions are clearly 
exponential [2] [69], others display noticeable curvature 
on a semi-log plot [47] [54] [70] [73] . Few simulations cap- 
ture more than three decades in the tail, making it dif- 
ficult to distinguish exponentials from broad Gaussians 
or compressed exponentials. The tail of P(f) may show 
signatures of the approach to isostaticity; this point is 
not settled, and the FNE, which vanishes at the isostatic 
point, provides no illumination. While data from Zhang 
and Makse cross over from Gaussian to exponential de- 
cay near unjamming [29], data from Silbert et al. remain 
Gaussian even when the distance to the critical packing 
fraction <j) c is as small as 10 -6 [73]. O'Hern et al. find 
Gaussian tails in an ensemble at fixed distance to the 
transition, while fluctuations in this distance due to finite 
size, which occur in fixed volume ensembles, can render 
tails exponential [70] . 
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FIG. 7: (a) Periodic force network on the frictionless trian- 
gular lattice, (b) Forces on each grain construct a tile (c) by 
summing them graphically. Tiles tessellate space due to New- 
ton's laws, (d) The reciprocal tiling of (a). (From Ref. [50].) 
(e) The numerical distribution P(p) of the pressure on a grain 
in the triangular lattice (solid curve) is well described by 
Eq. (|8| (black dashed curve) but not by a comparable cal- 
culation neglecting conservation of tiling area (gray dashed 
curve). (Adapted from [36].) 



III. MECHANICAL PROPERTIES 

It is generally thought that force networks are impor- 
tant for the mechanical properties of static and qua- 
sistatic granular materials. Since the FNE accurately 
describes the statistics of real force networks, one could 
ask whether it also captures mechanical response to ex- 
ternal loading. This can indeed be explored within the 
FNE. We discuss the response to anisotropic loads (shear 
stress) and localized loads (a point force at the bound- 
ary). 



A. Response to a shear stress 

Anisotropic force states such as shown in Fig. [5] appear 
naturally when imposing a shear stress to the system. It 
is interesting to study this effect for disordered packings, 
for which the anisotropy cannot be aligned along pre- 
ferred lattice directions. Figure |8^i illustrates how the 
imposed shear stress a xy gives rise to major and minor 
principal axes at angles <p = 7r/4 and 37r/4 respectively 
[42] . As above, we quantify anisotropy by r from Eq. [5] 

The resulting force anisotropy can be investigated us- 
ing the angle resolved PDF, representing P(f) at differ- 
ent orientations <p. Figure (8}d shows results obtained in 
the FNE for frictionless packings with isotropic contact 
networks [42] . The modulation along <j> increases with 
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FIG. 8: (a) Shear stress breaks isotropy, introducing a major 
and minor principal stress axis, (b) Density plot of the angle- 
dependent force distribution P</>(/) for increasing anisotropy 
t. (c) A strongly anisotropic force network. Adapted from 
Ref. 142]. 




FIG. 9: Photelastic granular packing under shear. The 
strongest forces align with the major principal axis. (From 
Ref. [75].) 



the imposed value of r, signaling the increase of force 
anisotropy. This is further quantified through the angle 
resolved force average, /(</>), which can be expanded in a 
Fourier series as in Refs. [42] 173]. 



f{<j>) = 1 + 2r sin 20 - b 2 cos 40 ■ 



(9) 



The second order coefficient is directly proportional to r, 
while the higher order terms do not couple to the stress 
at all. Indeed, Eq. ([9| already provides an accurate fit to 
the numerical results when truncated at fourth order. 

Force anisotropy has been observed experimentally [4, 
[75] and in contact dynamics simulations [74] . Figure [9] 
shows an experiment on sheared packings of photoela- 
sic grains, which visualizes the forces in the system [75] . 
These indeed reveal the preferred orientation of large 
forces along 45° with respect to the horizontal, coinciding 
with the major principal axis. 

Another striking feature of Fig. (8)3 is that many contact 
forces along the minor axis evolve towards zero for in- 
creasing shear (black area in Fig. 7b3-b4 near cf) ~ 37r/4). 
These small forces are close to breaking (they cannot be- 
come tensile) , which will eventually lead to failure of the 
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system. Thus there must be a maximum r m above which 
no solutions exist. One can interpret this in terms of the 
volume spanned by all force configurations, which contin- 
uously shrinks with r and reaches zero at r m . It is found 
numerically that r m strongly depends on the coordina- 
tion number 2, and approximately follows the scaling [42 

r m «2^^, (10) 

where Zi SO is the isostatic coordination number. One thus 
finds that the maximum possible stress vanishes when the 
packing approaches the isostatic limit. This is strongly 
reminiscent of the onset of bulk modulus, shear modulus 
and dynamic yield stress at the jamming transition. 

Far away from the isostatic point, it is possible to ap- 
ply mean field arguments to estimate the maximum shear 
stress in the FNE [76 . This is done by requiring the av- 
erage force /(</>) to be positive in all directions <p. This 
condition is of course much weaker than the requirement 
that all individual forces be positive, and therefore leads 
to an upper bound for r m . Generalizing the argument 
to frictional particles with friction coefficient /i, one esti- 
mates [76] 

< g • (11) 

Simulations of the FNE with z = 5.5 show that this upper 
bound is approached to within a few percent for /i = 0.5 
and 1.0 [47]. 

B. Response to a point force 

Another way to assess the mechanical behavior is 
through the response to a localized force. For small 
enough forces this probes the Green's function of the 
granular packing and provides crucial information on the 
effective continuum description of the system. Experi- 
ments [77-79 and simulations [3TJ [80j [81] have shown 
that the spreading of the load inside the material is not 
universal but can be along a single broad peak, as is 
the case for isotropic linear elastic materials, or more 
anisotropically along two peaks. The response is found 
to depend on parameters such as friction coefficient, de- 
gree of disorder, coordination number and amplitude of 
applied force. 

The response problem was addressed in the FNE for 
two-dimensional lattices with a free top surface [4lJ [43] . 
Frictionless grains were studied on a triangular lattice, 
while a square lattice was used for frictional systems. 
Here we discuss the findings of [43], where the packing 
was first put under isotropic pressure p before applying a 
load / on a grain in the top layer (Fig. [lO| . Within this 
setting it is possible to investigate the effect of the relative 
force amplitude, F = f/p, and the friction coefficient ja. 
The response function G(x, z) is defined by the difference 



in vertical force W(x,z) on a grain with and without the 
point force: 

G(x,z)= {W " )f - {W ** )o . (12) 
F 

The brackets with subscripts F and denote the average 
in the corresponding FNE. 

A linear regime was found for small enough load, 
F < 3, where the response is independent of the am- 
plitude F [43]. Figure [To| displays G(x,z) on the square 
lattice for various friction coefficients /i, all with F = 3. 
Each graph contains the response at different distances to 
the top surface and reveals how the load spreads through 
the material. Clearly, the response evolves from 'two 
peaks' to 'one peak' upon increasing \i. This can be inter- 
preted as follows. In the particular case where \i = 0, the 
response must be along the two downward lattice direc- 
tions. For fi 7^ 0, the presence of tangential forces makes 
it possible to spread the load. This 'freedom' increases 
with \i and eventually yields a single peaked response, as 
in isotropic elastic media. These observations agree well 
with MD simulations of ordered layers of grains [80] [82] , 
where the grain-grain interactions were modeled in detail. 
These simulations showed that friction strongly enhances 
the regime where the material behaves like a linear elastic 
solid. 

It is interesting to compare the \i — results for the 
square lattice (z = 4) and the triangular lattice (z = 6). 
Though it is dangerous to extrapolate generic properties 
of isostatic systems from lattice packings, it is notewor- 
thy that the response transitions from single-peaked (not 
shown) to double-peaked as (z — Zi SO ) — > + . This is 
consistent with the emerging picture that there is a di- 
verging length scale £* ~ 1/Az above which a packing 
can be viewed as an elastic continuum [67] . Anisotropy 
may also play a role. The triangular lattice is isotropic 
in linear elasticity while the square lattice is not, and 
anisotropic continua admit two-peaked response [83], [84] . 
Finally, at fixed depth the response gradually changes 
to two peaks when increasing F well into the nonlin- 
ear regime, also consistent with continuum descriptions 
[84] . This crossover occurs when the local load W be- 
comes much larger than the horizontal pressure scale, in 
which case the horizontal forces hardly contribute the 
force balance. Locally, this effectively changes the lat- 
tice from 'triangular' to 'square', the latter allowing only 
for transmission along the lattice directions - see also 
Ref. [80]. 

IV. CONCLUSIONS AND OUTLOOK 

The FNE is a convenient minimal model system for 
static granular materials that takes into account local 
force balance explicitly. Due to its simplicity its proper- 
ties can be computed accurately. 

The FNE reproduces the most robust features of the 
statistics of contact forces, suggesting that these features 
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FIG. 10: (a) Contact network and boundary conditions for 
point force response. An otherwise isotropic system is aug- 
mented by a normal force F at one boundary, (b) Evolution 
of the response with depth z changes qualitatively with mi- 
croscopic friction coefficient \i. (Adapted from 43 .) 

follow from the geometry of a high-dimensional space, 
and that details of the force law are of secondary impor- 
tance. Perhaps most surprising is the finding that P(f) 
in the FNE decays faster than exponentially. It is our im- 
pression that consensus has crystallized about the notion 
that exponential tails are a hallmark of granular force 
fluctuations, motivated by early experimental [lOl [TT] , 
numerical [2], and theoretical [12] work. Although there 
is support for this view, it seems to outrun the available 
evidence. The FNE suggests an alternative perspective; 
namely that distributions decay faster than exponentially 
asymptotically but, due to anisotropy, may appear expo- 
nential over accessible ranges. One useful role for the 
FNE is that of litmus test: a model that cannot explain 
results in the FNE is too simplistic. On this basis the 



q-model and ideal gas-like extrapolations from the Ed- 
wards ensemble can already be rejected. Therefore, if 
the tail of P(f) in real systems is robustly exponential, 
a theoretical explanation is still lacking. 

As a statistical measure, P(f) carries no information 
about the spatial structure of force networks. Studies 
of thresholded force networks offered the intriguing sug- 
gestion that systems with vector force balance represent 
a different universality class from ordinary percolation 
[44l [45]. More recent work, however, attributes these ob- 
servations to crossover effects [85] . It is therefore an open 
question whether force networks carry diverging spatial 
signatures of the impending loss of rigidity as the iso- 
static point is approached, as do the vibrational [66] and 
response [67] properties of soft sphere packings. 

Studies of response hint at a connection between hy- 
perstaticity of forces and elasticity or continuum-like re- 
sponse in the corresponding soft sphere packing [67]. It 
remains an interesting and open question whether the 
stress statistics of Eq. [8] break down near isostaticity. 

Important topics for future studies are the nature of 
force networks in frictional systems, in systems of non- 
spherical particles, and in flowing systems — does the 
FNE also capture the statistics of these systems? Fi- 
nally, one may wonder if the idea of flatly sampling a 
family of configurations, which the FNE does for force 
configurations, can be extended to other cases, such as 
the family of contact geometries that correspond to a 
certain contact topology. 



Acknowledgments 

BPT acknowledges support from the Dutch physics 
foundation FOM. 



[1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, 

Rev. Mod. Phys. 68, 1259 (1996). 
[2] F. Radjai, M. Jean, J. -J. Moreau, and S. Roux, 

Phys. Rev. Lett. 77, 274 (1996). 
[3] D. Howell, R. P. Behringer, and C. Veje, Phys. Rev. Lett. 

82, 5241 (1999). 
[4] T. S. Majmudar and R. P. Behringer, Nature 435, 1079 

(2005). 

[5] M. van Hecke, J. Phys. Cond. Matt, (accepted). 

[6] P. Dantu, in Proc. Of the J^th International Conf On 
Soil Mech. and Foundation Eng. (Butterworths Scientific 
Publications, London, 1957), vol. 1, pp. 144-148. 

[7] G. de Josselin de Jong and A. Verruijt, 
Cah. Gr. Franc. Rheol. 2, 73 (1969). 

[8] H. A. Janssen, Zeitschr. d. Vereines deutscher Ingenieure 
39, 1045 (1895). 

[9] M. Sperl, Gran. Matt. 8, 59 (2006). 
[10] C. H. Liu, S. R. Nagel, D. A. Schecter, S. N. Copper- 
smith, S. Majumdar, O. Narayan, and T. A. Witten, 
Science 269, 513 (1995). 



[11] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, 
Phys. Rev. E. 57, 3164 (1998). 

[12] S. N. Coppersmith, C. Liu, S. Majumdar, O. Narayan, 
and T. A. Witten, Phys. Rev. E 53, 4673 (1996). 

[13] L. Vanel, D. Howell, D. Clark, R. P. Behringer, and 
E. Clement, Phys. Rev. E 60, R5040 (1999). 

[14] M. Toiya, J. Stambaugh, and W. Losert, Phys. Rev. Lett. 
93, 088001 (2004). 

[15] R. M. Nedderman, Statics and Kinematics of Granu- 
lar Materials (Cambridge University Press, Cambridge, 
1992). 

[16] J. P. Wittmer, P. Claudin, M. Cates, and J.-P. Bouchaud, 

Nature 382, 336 (1996). 
[17] J. P. Bouchaud, in Slow Relaxations and Nonequilibrium 

Dynamics in Condensed Matter, edited by J. L. Barrat, 

M. Feigelman, J. Kurchan, and et al. (2004), pp. 131-197. 
[18] P. T. Metzger, Phys. Rev. E 77, 011307 (2008). 
[19] K. Brauer, M. Pfitzner, D. O. Krimer, M. Mayer, 

Y. Jiang, and M. Liu, Phys. Rev. E 74, 061311 (2006). 
[20] M. Depken, W. van Saarloos, and M. van Hecke, 



10 



Phys. Rev. E 73, 031302 (2006). 
[21] S. Henkes and B. Chakraborty, Phys. Rev. E 79, 061301 

(2009) . 

[22] J. H. Snoeijer, T. J. H. Vlugt, M. van Hecke, and W. van 
Saarloos, Phys. Rev. Lett. 92, 054302 (2004). 

[23] S. Alexander, Phys. Rep 296, 65 (1998). 

[24] C. F. Moukarzel, Phys. Rev. Lett. 81, 1634 (1998). 

[25] A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 
687 (1999). 

[26] J.-N. Roux, Phys. Rev. E 61, 6802 (2000). 

[27] A. J. Liu and S. R. Nagel, Nature 396, 21 (1998). 

[28] C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 

Phys. Rev. E 68, 011306 (2003). 
[29] H. P. Zhang and H. A. Makse, Phys. Rev. E 72, 011301 

(2005). 

[30] L. E. Silbert, G. S. Grest, and J. W. Landry, Phys. Rev. E 

66, 061303 (2002). 
[31] A. Kasahara and H. Nakanishi, Phys. Rev. E 70, 051309 

(2004). 

[32] T. Unger, J. Kertesz, and D. Wolf, Phys. Rev. Lett. 94, 
178001 (2005). 

[33] E. Somfai, M. van Hecke, W. G. Ellenbroek, 
K. Shundyak, and W. van Saarloos, Phys. Rev. E 75, 
060302 (R) (2007). 

[34] S. F. Edwards and R. B. S. Oakeshott, Physica A 157, 
1080 (1989). 

[35] J. H. Snoeijer, T. J. H. Vlugt, W. G. Ellenbroek, M. van 
Hecke, and J. M. J. van Leeuwen, Phys. Rev. E 70, 
061306 (2004). 

[36] B. P. Tighe, A. R. T. van Eerd, and T. J. H. Vlugt, 

Phys. Rev. Lett. 100, 238001 (2008). 
[37] S. F. Edwards, J. Phys. A 41, 324019 (2008). 
[38] T. C. Halsey and D. Ertas, Phys. Rev. Lett 83, 5007 

(1999). 

[39] S. McNamara and H. Herrmann, Phys. Rev. E 70, 061303 

(2004) . 

[40] B. P. Tighe, J. E. S. Socolar, D. G. Schaeffer, W. G. 
Mitchener, and M. L. Huber, Phys. Rev. E 72, 031306 

(2005) . 

[41] S. Ostojic and D. Panja, Europhys. Lett. 71, 70 (2005). 
[42] J. H. Snoeijer, W. G. Ellenbroek, T. J. H. Vlugt, and 

M. van Hecke, Phys. Rev. Lett. 96, 098001 (2006). 
[43] S. Ostojic and D. Panja, Phys. Rev. Lett. 97, 208001 

(2006) . 

[44] S. Ostojic, E. Somfai, and B. Nienhuis, Nature 439, 828 

(2006) . 

[45] S. Ostojic, T. J. H. Vlugt, and B. Nienhuis, Phys. Rev. E 

75, 030301(R) (2007). 
[46] A. R. T. van Eerd, W. G. Ellenbroek, M. van Hecke, J. H. 

Snoeijer, and T. J. H. Vlugt, Phys. Rev. E 75, 060302(R) 

(2007) . 

[47] A. R. T. van Eerd, B. P. Tighe, and T. J. H. Vlugt, 

Molecular Simulation 35, 1029 (2009). 
[48] M. R. Shaebani, T. Unger, and J. Kertesz, Phys. Rev. E 

79, 052302 (2009). 
[49] B. P. Tighe, in Powders and Grains 2009, edited by 

M. Nakagawa and S. Luding (American Institute of 

Physics, 2009), pp. 305-308. 
[50] B. P. Tighe and T. J. H. Vlugt, J. Stat. Mech. P01015 

(2010) . 

[51] C. S. O'Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, 

Phys. Rev. Lett. 86, 111 (2001). 
[52] H. A. Makse, D. L. Johnson, and L. M. Schwartz, 

Phys. Rev. Lett. 84, 4160 (2000). 



[53] A. V. Tkachenko and T. A. Witten, Phys. Rev. E 62, 
2510 (2000). 

[54] C. Goldenberg and I. Goldhirsch, Granular Matter 6, 87 

(2004) . 

[55] O. Tsoungui, D. Vallet, and J.-C. Charmet, Granular 

Matter 1, 65 (1998). 
[56] G. L0voll, K. J. Mal0y, and E. G. Flekk0y, Phys. Rev. E 

60, 5872 (1999). 
[57] J. Brujic, S. F. Edwards, I. Hopkinson, and H. A. Makse, 

Physica A 327, 201 (2003). 
[58] J. Zhou, S. Long, Q. Wang, and A. D. Dinsmore, Science 

312, 1631 (2006). 
[59] J. C. Maxwell, Philosoph. Mag. 27, 250 (1864). 
[60] P. Evesque, Poudres et grains 9, 13 (1999). 
[61] N. P. Kruyt and L. Rothenburg, Int. J. Solids Structures 

39, 571 (2002). 
[62] K. Bagi, Granular Matter 5, 45 (2003). 
[63] A. H. W. Ngan, Phys. Rev. E 68, 011301 (2003). 
[64] J. D. Goddard, Int. J. Solids Structures 41, 5851 (2004). 
[65] S. Henkes, C. S. O'Hern, and B. Chakraborty, 

Phys. Rev. Lett. 99, 038002 (2007). 
[66] M. Wyart, S. R. Nagel, and T. A. Witten, Euro- 
phys. Lett. 72, 486 (2005). 
[67] W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van 

Saarloos, Phys. Rev. Lett. 97, 258001 (2006). 
[68] R. Schneider, in Handbook of convex geometry, edited by 

P. Gruber and J. M. Wills (North Holland, Amsterdam, 

1993), pp. 273-299. 
[69] F. Radjai, S. Roux, and J.-J. Moreau, Chaos 9, 544 

(1999). 

[70] C. S. O'Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, 

Phys. Rev. Lett. 88, 075507 (2002). 
[71] D. L. Blair, N. W. Mueggenburg, A. H. Marshall, H. M. 

Jaeger, and S. R. Nagel, Phys. Rev. E 63, 041204 (2001). 
[72] J. M. Erikson, N. W. Mueggenburg, H. M. Jaeger, and 

S. R. Nagel, Phys. Rev. E 66, 040301(R) (2002). 
[73] L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E p. 

041304 (2006). 

[74] F. Radjai, D. E. Wolf, M. Jean, and J. J. Moreau, 

Phys. Rev. Lett. 80, 61 (1998). 
[75] A. P. F. Atman, P. Brunet, J. Geng, G. Reydel- 

let, P. Claudin, R. P. Behringer, and E. Clment, 

Eur. Phys. J. E 17, 93 (2005). 
[76] W. G. Ellenbroek and J. H. Snoeijer, J. Stat. Mech. p. 

P01023 (2007). 

[77] J. Geng, D. Howell, E. Longhi, R. P. Behringer, 

G. Reydellet, L. Vanel, E. Clement, and S. Luding, 

Phys. Rev. Lett. 87, 035506 (2001). 
[78] D. Serero, G. Reydellet, P. Claudin, E. Clement, and 

D. Levine, Eur. Phys. J. E 6, 169 (2001). 
[79] N. W. Mueggenburg, H. M. Jaeger, and S. R. Nagel, 

Phys. Rev. E 66, 031304 (2002). 
[80] C. Goldenberg and I. Goldhirsch, Phys. Rev. Lett. 89, 

084302 (2002). 

[81] N. Gland, P. Wang, and H. A. Makse, Eur. Phys. J. E 
20, 179 (2006). 

[82] C. Goldenberg and I. Goldhirsch, Nature 435, 188 

(2005) . 

[83] M. Otto, J.-P. Bouchaud, P. Claudin, and J. E. S. Soco- 
lar, Phys. Rev. E 67, 031302 (2003). 

[84] B. P. Tighe and J. E. S. Socolar, Phys. Rev. E 77, 031303 
(2008). 

[85] B. Nienhuis (private communication). 



